




“The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text.” – jupyter.org

import sys
v = sys.version_info
print(f"We are using Python {v.major}.{v.minor}.{v.micro}")
We are using Python 3.9.5
The output can also be interactive.
from ipywidgets import interact
def f(x):
return x
interact(f, x=10);
interactive(children=(IntSlider(value=10, description='x', max=30, min=-10), Output()), _dom_classes=('widget-…
GIS714: geospatial computation and simulation: graduate level course required for Geospatial Analytics PhD students

You can try it out in Binder! (also linked from the GRASS GIS GitHub README)
# Import Python standard library and IPython packages we need.
import subprocess
import sys
# Ask GRASS GIS where its Python packages are.
sys.path.append(
subprocess.check_output(["grass80", "--config", "python_path"], text=True, shell=True).strip()
)
# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj
# Start GRASS Session
session = gj.init("data", "nc_spm_08_grass7", "PERMANENT")
# Set region to the high resolution study area
gs.run_command("g.region", region="rural_1m")
And then we are ready to begin using GRASS!
We can use the command line interface for GRASS by using the cell magic ! (for individual lines) or %%bash (applies to entire cell) which sends the code to the terminal.
!g.list type=raster -m -t
raster/aspect@PERMANENT raster/basin_50K@PERMANENT raster/boundary_county_500m@PERMANENT raster/cfactorbare_1m@PERMANENT raster/cfactorgrow_1m@PERMANENT raster/el_D782_6m@PERMANENT raster/el_D783_6m@PERMANENT raster/el_D792_6m@PERMANENT raster/el_D793_6m@PERMANENT raster/elev_lid792_1m@PERMANENT raster/elev_ned_30m@PERMANENT raster/elev_srtm_30m@PERMANENT raster/elev_state_500m@PERMANENT raster/elevation@PERMANENT raster/elevation_shade@PERMANENT raster/facility@PERMANENT raster/geology_30m@PERMANENT raster/lakes@PERMANENT raster/landclass96@PERMANENT raster/landcover_1m@PERMANENT raster/landuse96_28m@PERMANENT raster/lsat7_2002_10@PERMANENT raster/lsat7_2002_20@PERMANENT raster/lsat7_2002_30@PERMANENT raster/lsat7_2002_40@PERMANENT raster/lsat7_2002_50@PERMANENT raster/lsat7_2002_61@PERMANENT raster/lsat7_2002_62@PERMANENT raster/lsat7_2002_70@PERMANENT raster/lsat7_2002_80@PERMANENT raster/ncmask_500m@PERMANENT raster/ortho_2001_t792_1m@PERMANENT raster/roadsmajor@PERMANENT raster/slope@PERMANENT raster/soilsID@PERMANENT raster/soils_Kfactor@PERMANENT raster/streams_derived@PERMANENT raster/towns@PERMANENT raster/urban@PERMANENT raster/zipcodes@PERMANENT raster/zipcodes_dbl@PERMANENT
Or we can use the Python API for GRASS GIS. Since grass.jupyter is written in Python, we'll use the GRASS Python API from here on out.
gs.read_command("g.list", type="vector", flags="mt")
'vector/P079214@PERMANENT\r\nvector/P079215@PERMANENT\r\nvector/P079218@PERMANENT\r\nvector/P079219@PERMANENT\r\nvector/boundary_county@PERMANENT\r\nvector/boundary_municp@PERMANENT\r\nvector/bridges@PERMANENT\r\nvector/busroute1@PERMANENT\r\nvector/busroute11@PERMANENT\r\nvector/busroute6@PERMANENT\r\nvector/busroute_a@PERMANENT\r\nvector/busroutesall@PERMANENT\r\nvector/busstopsall@PERMANENT\r\nvector/census_wake2000@PERMANENT\r\nvector/censusblk_swwake@PERMANENT\r\nvector/comm_colleges@PERMANENT\r\nvector/elev_lid792_bepts@PERMANENT\r\nvector/elev_lid792_cont1m@PERMANENT\r\nvector/elev_lid792_randpts@PERMANENT\r\nvector/elev_lidrural_mrpts@PERMANENT\r\nvector/elev_lidrural_mrptsft@PERMANENT\r\nvector/elev_ned10m_cont10m@PERMANENT\r\nvector/firestations@PERMANENT\r\nvector/geodetic_pts@PERMANENT\r\nvector/geodetic_swwake_pts@PERMANENT\r\nvector/geology@PERMANENT\r\nvector/geonames_NC@PERMANENT\r\nvector/geonames_wake@PERMANENT\r\nvector/hospitals@PERMANENT\r\nvector/lakes@PERMANENT\r\nvector/nc_state@PERMANENT\r\nvector/overpasses@PERMANENT\r\nvector/poi_names_wake@PERMANENT\r\nvector/precip_30ynormals@PERMANENT\r\nvector/precip_30ynormals_3d@PERMANENT\r\nvector/railroads@PERMANENT\r\nvector/roadsmajor@PERMANENT\r\nvector/schools_wake@PERMANENT\r\nvector/soils_general@PERMANENT\r\nvector/soils_wake@PERMANENT\r\nvector/streams@PERMANENT\r\nvector/streets_wake@PERMANENT\r\nvector/swwake_10m@PERMANENT\r\nvector/urbanarea@PERMANENT\r\nvector/usgsgages@PERMANENT\r\nvector/zipcodes_wake@PERMANENT\r\n'
Add layers:
d.rast -> map.d_rast(), map3D.d_rast
show and save
folium...

Moving data from GRASS GIS location projection to folium is a challenge
folium is projected in Web Mercator (EPSG:3857)
However, any coordinates (i.e. vector data) need to be specified in degrees of latitude and longitude (WGS84, EPSG:4326)
folium reprojects latitude and longitude to Web Mercator internally
We pass data to folium by reprojecting to temporary locations

import folium
# Create figure
fig = folium.Figure(width=600, height=400)
# Create a map to add to the figure later
m = folium.Map(tiles="Stamen Terrain", location=[35.761168,-78.668271], zoom_start=13)
# Create and add elevation layer to map
gj.Raster("elevation", opacity=0.5).add_to(m)
# Do some cool folium stuff!
# Like make a tooltip
tooltip = "Click me!"
# and add a marker
folium.Marker(
[35.781608,-78.675800], popup="<i>Center For Geospatial Analytics</i>", tooltip=tooltip
).add_to(m)
# and a circle
folium.Circle(
radius=120,
location=[35.769781,-78.663160],
popup="Great Picnic Area",
color="crimson",
fill=False,
tooltip=tooltip
).add_to(m)
# Add the map to the figure
fig.add_child(m)
# Display figure
fig
The InteractiveMap class allows users unfamiliar with folium to produce maps easily.
# Create Interactive Map
fig = gj.InteractiveMap(width = 600)
# Add raster, vector and layer control to map
fig.add_raster("elevation", title="Elevation Raster", opacity=0.8)
#fig.add_vector("roadsmajor")
fig.add_layer_control(position = "bottomright")
# Display map
fig.show()
GIS714: Geospatial Computation and Simulation, spring 2022,
# Make a new mapset for this assignment
gs.run_command("g.mapset", mapset="HW3_water_simulation", location="nc_spm_08_grass7", flags="c")
We create a map of flooded area with r.lake (GRASS manual for r.lake) by providing a water level and a seed point:
gs.run_command("r.lake", elevation="elev_lid792_1m", water_level=113.5, lake="flood1", coordinates="638728,220278")
# See results
flood1 = gj.Map(use_region=True)
flood1.d_rast(map="elev_lid792_1m")
flood1.d_rast(map="flood1")
# Display map
flood1.show()
Increase water level to 113.7m and 114.0m and create flooded area maps at these two levels.
#### Your Answer Here
We will use two GRASS addons, r.stream.distance and r.lake.series, to estimate inundation with Height Above Nearest Drainage methodology (A.D. Nobre, 2011). We need to install them first:
#!g.extension r.stream.distance
#!g.extension r.lake.series
For this section, we will change our computation region to elevation which is a larger study area than we used above. Because we are using a new region (and also a higher threshold of 100,000), we need to run r.watershed again to compute the flow accumulation, drainage and streams. We convert the streams to vector for better visualization.
gs.run_command("g.region", raster="elevation")
gs.run_command("r.watershed", elevation="elevation", accumulation="flowacc", drainage="drainage", stream="streams_100k", threshold=100000)
gs.run_command("r.to.vect", input="streams_100k", output="streams_100k", type="line")
Now we use r.stream.distance with output parameter difference to compute new raster where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains. This is our HAND terrain model.
gs.run_command("r.stream.distance", stream_rast="streams_100k", direction="drainage", elevation="elevation", method="downstream", difference="above_stream")
With r.lake.series, we can create a series of inundation maps with rising water levels. r.lake.series creates a space-time dataset. We can use temporal modules to further work with the data...
gs.run_command("r.lake.series", elevation="above_stream", start_water_level=0, end_water_level=5, water_level_step=0.5,
output="inundation", seed_raster="streams")
... or visualize the flood with TimeSeriesMap.
flood_map = gj.TimeSeriesMap(use_region=True)
flood_map.add_raster_series("inundation")
flood_map.d_legend(color="black", at=(10,40,2,6)) #Add legend
flood_map.d_barscale()
flood.show()
Notebook format was generally well received:
There were several areas of the notebook deployment that needed improvement or that posed a difficulty to the class:
FOSS4G '22 GRASS GIS Workshop Materials (Anna Petrasova)
GIS714: Advanced Geospatial Computation and Simulation Course Repository
Helena Mitasova
Stefan Blumentrath
Vero Andreo
GIS714 Class Spring 2021
… and the GRASS community, NC State Center for Geospatial Analytics, Google Summer of Code